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I.  ACCOMPLISHMENTS 

The  following  is  a  brief  summary  of  the  research  accomplishments  funded  by  AFOSR-88-0262,  a  grant  to 
seven  researchers  of  the  Institute  for  Computatonal  Mathematics  and  Applications,  University  of  Pittsburgh. 

1.  Additive  Correction  Methods. 

Using  a  general  theorem  that  gives  a  bound  on  the  iteration  error  in  terms  of  an  approximation  error  for 
the  subspace  containing  the  additive  correction,  realistic  mesh  independent  contraction  numbers  have  been 
derived  for  the  convergence  of  a  model  two  grid  method.  A  second  result  shows  that  if  the  subspaces  are  chosen 
to  be  certain  Krylov  subspaces  related  to  pseudo-residual  vectors  produced  by  the  iterations,  then  the  method  is 
exactly  a  restarted  preconditioned  conjugate  gradient  method  in  which  the  splitting  matrix  plays  the  role  of  the 
preconditioner.  Thus,  the  additive  correction  methods  include  this  important  class  of  iterative  methods  as  a  spe¬ 
cial  case. 

Regarding  the  convergence  of  generic  iterative  methods  employing  an  additive  correction  phase,  it  has 
been  shown  [1]  that  if  the  coefficient  matrix  is  symmetric,  positive  definite  (SPD),  if  each  additive  correction 
step  uses  orthogonal  projection  with  respect  to  the  "energy"  inner  product,  and  if  the  range  of  each  projector 
contains  the  most  recent  residual  vector,  then  a  contraction  number  for  the  error  is  (k  -  1)/(k  +  1),  where  k  is 
the  spectral  condition  number  of  the  system.  Thus,  such  methods  are  necessarily  convergent.  This  work  has 
recently  been  generalized  [2]  to  the  case  where  the  coei'ficU-t  matrix  is  only  positive  real;  i.e.,  its  symmetric 
part  is  SPD.  The  analysis  in  [2]  is  also  general  enough  to  yield  contraction  numbers  for  restarted  versions  of  the 
ORTHOMIN  and  GCW  generalized  conjugate  gradient  methods  described  in  [3]. 

It  appears  possible  to  sharpen  the  contracuon  numbers  obtained  in  [1]  and  [2]  by  incorporating  more  infor¬ 
mation  on  the  range  of  the  projectors.  This  leads  to  the  notion  of  abstract  angles  in  projection  methods. 

References 

[1]  Chou,  S.  H.  and  Porsching,  T.  A.,  "A  Note  on  Contraction  Numbers  for  Additive  Correction  Methods", 
Appl.  Math.  Lett.  2  (1989),  83-86. 

[2]  Chou,  S.  H.  and  Porsching,  T.  A.,  "Contraction  Numbers  for  Additive  Correction  Methods",  to  appear  Lin. 
Alg.  and  its  Appl.. 

13]  Hageman,  L.  A.  and  Young,  D.  M.,  Applied  Iterative  Methods,  Academic  Press,  New  York,  1981. 
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2.  Iterative  Solution  to  Navier-Stokes  Difference  Equations. 

The  time-dependent  two-dimensional  Navier-Stokes  equations  are  used  to  model  the  evolution  of  the  flow 
of  a  Newtonian  fluid.  Implicit  finite  difference  equations  on  a  MAC  discretization  grid  are  used  to  approximate 
the  continuity,  enthalpy  and  momentum  equations.  By  time-lagging  the  pressure  and  velocity  variables  in  the 
enthalpy  equations,  the  discrete  enthalpy  equations  are  uncoupled  from  the  continuity  and  momentum  equations 
and  solved  separately.  The  coupled  system  of  discrete  momentum  and  continuity  equations  is  then  transformed 
into  the  dual  variable  system  [1,  2,  3],  which  is  one-third  the  size  of  the  coupled  system. 

The  dual  variable  system  was  shown  to  be  solvable  for  all  problems  unless  the  time  step  chosen  is  one  of 
a  finite  number  of  values.  New  conditions  which  guarantee  the  solvability  of  the  system  for  all  positive  real 
values  of  the  time  step  were  derived. 

Special  iterative  methods  for  solving  the  dual  variable  system  were  developed.  These  methods  were 
shown  to  converge  to  the  solution  under  a  variety  of  conditions.  Also  techniques  for  accelerating  these  iterative 
methods  were  constructed.  These  special  iterative  methods  require  the  solution,  on  each  iteration,  of  a  linear 
system  whose  structure  resembles  a  discrete  Laplace  equation;  special  methods  of  handling  this  with  a  minimum 
of  operations  were  investigated. 

DUALIT  is  a  computer  code  which  solves  the  time-dependent  two-dimensional  Navier-Stokes  problem.  To 
solve  the  dual  variable  system,  DUALIT  uses  the  transformed  Jacobi  iterative  method  and  accelerates  conver¬ 
gence  with  a  second-order  Richardson  method.  The  aforementioned  discrete  Laplace  system  is  solved  with 
reduced-system  Richardson  iteration  accelerated  by  the  conjugate  gradient  method  [4],  A  comparison  of 
DUALIT  with  DUVAL,  a  commercial  program  which  uses  banded  Gaussian  elimination  to  solve  the  same  sys¬ 
tem  of  equations,  shows  that  the  dual  variable  solver  in  DUALIT  uses  less  computer  memory  and  considerably 
less  computer  time  to  solve  the  same  problems.  In  [7],  a  simulation  of  the  flow  of  air  along  the  exterior  of  an 
aircraft  and  into  an  opening  in  the  fuselage  is  anlayzed  numerically.  Inside  the  cavity,  there  is  a  sensor  which  is 
treated  as  a  blockage.  In  front  of  the  opening  in  the  fuselage  is  a  ramp  or  spoiler  which  causes  the  jetstream  of 
air  to  shoot  up  over  the  opening  and  reduces  the  amount  of  flow  into  the  cavity. 


The  inlet  velocities  on  the  bottom  of  the  region  correspond  to  a  free  stream  Mach  number  of  .75  while 


ttt 


Figure  i.  The  Aircraft  Cavity 

pressures  of  14.7  psi  are  specified  at  the  outlet  on  the  top  of  the  region.  The  walls  of  the  cavity  and  sensor  are 
no-slip  walls.  The  DUALIT  code  cannot  handle  walls  on  the  interior  of  the  region;  therefore,  the  fuselage, 
spoiler,  and  sensor  are  approximated  by  adding  large  friction  factors  to  the  appropriate  momentum  equations  to 
force  the  mass-velocities  to  be  zero. 

The  solutions  achieved  by  DUALIT  and  DUVAL  agree  to  2  significant  places  throughout  the  transient. 
Velocity  magnitude  contours  are  given  in  Figure  2  at  time  25  seconds.  From  the  figure  one  .in  see  the 
effectiveness  of  using  friction  factors  to  simulate  the  interior  blockages. 

Figure  3  shows  the  streamlines  of  the  flow  at  25.0  seconds.  As  expected,  the  spoiler  causes  most  of  the 
air  to  stream  up  and  over  the  opening.  There  is  a  small  vortex  just  above  and  to  the  left  of  the  sensor. 

Figure  4  is  representative  of  the  gain  in  efficiency  using  the  algorithm  in  DUALIT. 

This  research  is  contained  in  a  Fh.D.  thesis  [5]  of  George  Mesina,  a  student  of  Professor  Charles  Hall,  and 


is  also  contained  in  a  forthcoming  joint  paper  [6], 


Figure  2.  Aircraft  Cavity  Velocity  Magnitude  Contours  at  25.0  Seconds 


We  have  developed  an  algorithm  to  construct  cycle  bases  for  MAC  finite  difference  equations  to  extend 
the  dual  variable  method  [1,  2]  to  three  dimensional  analyses.  This  has  been  implemented  on  the  CRAY/XMP. 
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Figure  4.  Comparison  of  Iterative  and  Direct  Solvers  on  the  CRAY/XMP48 
An  initial  testing  phase  has  been  completed  in  which  uniform  flow,  PoiseuiUe  flow  and  driven  cavity  flows  have 
been  simulated  successfully. 
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3.  Solving  Discrete  Navier-Stokes  Systems. 

The  focus  of  this  project  was  solving  the  system  of  equations  arising  from  discretizations  of  the 
incompressible  Navier-Stokes  equations  in  an  efficient  way.  This  remains  one  of  the  most  challenging  problems 
in  computational  fluid  dynamics,  with  at  least  thr’e  basic  issues: 

(a)  Linearization  of  the  nonlinear  system, 

(b)  Handling  the  incompressibility  constraint, 

(c)  Finding  an  alternative  less  costly  way  to  solve  convection  dominated  diffusion  equations  than  direct 
methods. 

Our  approach  to  dealing  with  these  issues  was  to  use  a  new  iterative  method  based  on  conjugate  gradient 
(eg  )  techniques.  This  was  done  with  the  aim  of  preserving  the  attractive  features  of  eg ,  namely: 

(i)  Small  storage  requirements:  O  (N)  versus  0  (A53)  for  efficient  banded  storage  for  direct  methods  needed 
to  solve  2nd  order  partial  differential  equations  in  3  space  dimensions. 

(ii)  Low  computational  cost:  0(V4'3)  versus  0(/V7,3)  for  efficient  banded  elimination. 

(iii)  Ease  of  implementation:  No  parameters  need  to  be  supplied  by  the  user,  as  for  S.O.R.  Also,  it  is  not  lim¬ 
ited  to  simple  geometry,  as  is  multi-grid. 

(iv)  Convergence:  eg  iterates  are  theoretically  proven  to  their  error  decreasing  monotonically,  and  conver¬ 
gence  is  guaranteed  within  N  iterations.  (In  practice,  the  iterations  may  be  stopped  on  before  N ,  since 
sufficient  accuracy  is  obtained  after  only  O  (V 1/3 log  Nlt3)  iterations. 

Unfortunately,  the  eg  method  was  designed  for  linear  symmetric  systems  of  equations  Mx  =  b.  Thus  it 
could  not  be  applied  to  Navier  Stokes  disci aizuiions,  which  are  nonsymmetric  nonlinear  systems  with  additional 


(incompressibility)  constraints.  Bi-conjugate  gradients  (beg )  is  a  generalization,  designed  for  unsymmetric 
matrices.  In  describing  how  beg  arises  from  eg,  we  need  to  use  a  three-term  recurrence  relation  satisfied  by  eg 
iterates  xk 

*k+i  =  <£kxk  +  (1  -  tOtjr*. |  +  to kctkrk  . 

This  is  not  the  more  standard  description  of  eg,  which  is  in  terms  of  conjugate  directions,  as  in  [2].  It  follows 
that  the  residuals  rk  =  b  -  Mxk  also  satisfy  a  similar  relation: 

rk+ 1  =  <-■->* rk  +  (1  -  to* )rk-\  ~  U>ka.kMrk  . 

Beg  introduces  a  new  sequence  of  %'ectors  defined  tu  satisfy  the  residua!  relation  for  the  transpose  of  the  matrix: 

Sk+ i  +  (1  -  1  . 

The  coefficients  a*  and  w*  are  defined  as  follows 

_  <St  k> 
a*  <5* ,  Mrk  > 

I  ^  «t  <Sk,  rk>  1 
COt  o.k-\  <Sk.\,rk.i^  u,.| 

Beg  preserves  the  (i)  storage,  (ii)  computational  cost,  and  (iii)  algorithmic  advantages  of  eg ,  but  loses  the 
theoretical  convergence  property  (iv).  In  fact,  the  residual  norms  fluctuate  quite  drasticaly,  making  the  method 
unusable.  To  remedy  this,  we  introduce  smoothing  via  an  auxiliary  sequence  pt  of  averaged  residuals,  which 
correspond  to  averaged  solutions  c* . 

The  smoothing  forces  monotonic  descent  of  residual  norms,  making  beg  a  usable  method.  The  cost  of 
implementing  smoothing  is  ncgligable,  the  additional  storage  being  4  vectors,  and  the  additional  computations 
being  2  scalar  multiplications  of  vectors  and  ?  dot  products.  Furthermore,  the  smoothing  is  only  a  "post¬ 
processing"  operation,  without  interfering  in  any  way  with  the  beg  algorithm  itself. 

Brg  with  smoothing  is  a  computationally  cheaper  iterative  method  for  nonsymmctric  linear  systems  than 
others  such  as  Orthonium  [1],  both  in  terms  of  op-counts  and  storage  required.  It  constitutes  one  approach  to 
issue  (c)  in  the  context  of  Navicr-Stokcs,  We  deal  with  issue  (b)  by  actualy  carrying  out  the  beg  algorithm  in 
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the  space  of  solenoidal  vector  fields  bv  means  of  projections.  Issue  (a)  is  the  easiest  to  handle,  via  lagging. 
Substantial  numerical  computations  with  this  new  algorithm  are  reported  in  [3], 

References 

[1]  H.  Elman,  Yale  Univ.  Report. 

[2]  D.  G.  Lucnbergcr  (1984):  Linear  and  Nonlinear  Programming,  2nd  edition,  Addison-Wesley,  p.  244. 

[3]  S.  Choudhury,  "A  Projected  Bi-Conjugate  Gradient  Approach  to  Solving  Discretized  Incompresible 
Navier-Stokes  Equations",  ICMA  Report. 

4.  Divergence  Free  Subspaces. 

Finite  element  algorithms  for  solving  the  incompressible  fluid  flow  problems  must  deal  with  the  conserva¬ 
tion  of  mass  equations: 

div  g  =  0  (1) 

where  q  is  the  mass  velocity  vector.  Penality  methods  seem  to  be  the  popular  approach,  although  other  investi¬ 
gators  have  been  able  to  construct  special  divergence  free  basis  elements  that  satisfy  (1)  a  priori. 

Knowing  such  a  basis,  Galcrkin’s  method  involves  a  projection  into  this  reduced  subspace  of  divergence 
free  velocity  fields.  Unfortunately,  efforts  in  this  direction  up  to  now  have  involved  straight  sided  elements  and 
low  order  element  approximation  to  the  pressure. 

The  dual  variable  method,  first  introduced  by  Amil,  Hall  and  Porsching  [1],  finds  a  basis  for  the  null  space 
of  the  discrete  divergence  operator  using  network  theory.  It  was  applied  to  the  finite  difference  scheme  pro¬ 
posed  by  Kr/.hivitski  and  Ladyzhenskaya  (2]  and  a  scheme  using  the  discrete  divergence  and  gradient  operators 
of  the  MAC  method  studied  by  Harlow  and  Wcich  [3].  The  DVM  decouples  the  pressure  from  the  velocity  and 
results  in  a  large  reduction  in  the  size  of  the  linear  system  to  be  solved  at  each  time  step.  Dougall,  Hall  and 
Porsching  [4 1  implemented  this  method  for  finite  difference  discretization  of  Navier-Stokes  equations  used  in  the 
simulation  of  (low  of  thermally  expandable  steam  water  mixtures  in  reactor  components. 

Cha  and  Porsching  [51  extended  this  method  to  finite  difference  discretizations  of  steady-state  Navier- 
Stokes  equations.  Sledge  [61  and  Hall,  Peterson,  Porsching  and  Sledge  [7]  applied  the  DVM  to  finite  element 
discretizations  of  transient  Navier-Stokes  equations.  This  work  involved  a  4-node  quadrilateral  clement  for 
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velocity  and  a  constant  quadrilateral  element  for  pressure.  No  elements  with  curved  edges  were  considered. 
Mansulti,  et  al  [8]  and  Bulgarelli,  et  al  [9]  applied  the  dual  variable  method  to  two  and  three  dimensional  prob¬ 
lems  using  finite  difference  approximations.  Burkardt,  Hall  and  Porsching  [101  formulate  a  DVM  for  two 
dimensional  compressible  flow  problems  and  Frey,  Hall  and  Porsching  [11]  apply  the  DVM  to  analyze  confined 
aerodynamical  flows.  Goodrich  and  Soh  [12]  interpreted  the  null  basis  in  terms  of  the  stream  function  for  2-D 
incompressible  flow  in  driven  cavities. 

There  have  been  many  constructions  of  explicit  bases  of  the  divergence  free  subspace  for  various  finite 
element  and  finite  difference  schemes.  Griffiths  [13,  14,  15)  obtained  an  clement  level  divergence  free  basis  for 
several  finite  element  schemes  on  triangular  and  quadrilateral  elements.  The  techniques  used  are  the  same  in 
each  paper.  Approximate  values  of  the  stream  function  at  comer  nodes  are  used  to  eliminate  the  unknown  velo¬ 
city  components  at  midside  nodes  so  that  a  :_.pi,.ul  divergence  free  function  on  each  element  is  derived. 

In  Griffiths  [13,  14],  three  types  of  finite  element  schemes  were  investigated  on  triangular  elements  which 
were  given  by  Crouzeix  and  Raviart  [16],  A  divergence  free  basis  was  given  for  a  nonconforming  velocity  field 
where  the  components  of  velocity  are  represented  by  piccewfisc  linear  functions  defined  in  terms  of  their  values 
at  the  midside  nodes  of  the  triangles.  A  divergence  free  basis  also  was  given  lor  a  velocity  field  where  the  com¬ 
ponents  of  the  velocity  are  piecewise  quadratic  functions  defined  in  terms  of  values  at  the  vertices  and  midside 
nodes  of  each  triangle.  The  (discontinuous)  piecewise  constant  pressure  space  was  used  for  both  of  the  above 
velocity  spaces.  Another  divergence  free  subspace  derived  by  Griffiths  [13]  involved  a  velocity  field  which 
comes  from  adding  a  cubic  term  to  the  quadratic  representation.  The  pressure  space  used  was  a  piecewise  linear 
function  with  a  single  element  support  Griffiths  [15]  derived  a  basis  for  the  divergence  free  subspace  of  the  9- 
node  biquadratic  element  velocity  field  on  quadrilatcrial  elements.  The  following  corresponding  pressure  spaces 
were  investigated:  constant,  linear  and  bilinear  elements.  The  basis  functions  for  these  pressure  spaces  have  sup¬ 
port  on  a  single  element.  This  allows  the  incompressible  constraint  to  be  analyzed  one  clement  at  a  time.  But 
the  basis  function  for  the  pressure  space  is  discontinuous  at  clement  boundaries. 

Gustafson,  et  al  [17]  combine  group  theory  and  fluid  mechanics  theory  to  obtain  a  basis  for  the  divergence 
free  subspace  associated  with  the  choice  of  quadratic  velocity  and  constant  pressure  triangular  elements  in  two 
dimension.  Similar  results  have  been  obtained  in  3-D  for  the  scheme  referred  to  as  APX  3  in  Teman  [18].  The 
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late  work  by  Gustufson,  et  al  [19]  can  be  viewed  as  augmenting  and  ending  their  previous  work. 

Stephens,  et  al  [20]  and  Goodrich  and  Soh  [12]  applied  the  Galerkin  finite  difference  method  (GFDM)  to 
the  incompressible  Navier-Stokes  problem.  This  approach  is  similar  to  the  Galerkin  finite  element  method.  The 
subspaces  of  discrete  divergence  free  mesh  vectors  were  constructed  for  several  finite  difference  schemes.  It  is 
required  that  the  discrete  divergence  and  discrete  gradient  operators  are  formally  adjoint.  In  Fortin  [21],  a  sub¬ 
space  of  Sh  is  con  meted  in  which  a  function  satisfies  (1)  for  a  subspace  of  the  pressure  space.  This  subspace 
is  the  orthogonal  complement  of  piecewise  constant  pressure  space.  This  can  reduce  the  5-node  velocity  and 
linear  discontinuous  pressure  element  to  4-node  velocity  and  discontinuous  constant  pressure  element. 

We  emphasize  that  all  of  the  above  constructions  of  divergence  free  basis  velocity  vectors  require  straight 
sided  elements  and  that  the  pressure  space  contains  only  dicontinuous  functions. 

For  many  finite  element  and  finite  difference  schemes,  divergence  free  bases  and  null  bases  have  been 
obtained.  But  most  of  them  are  for  low  order  schemes  and  for  very  regular  domains  where  the  boundary  is  a 
piecewise  straight  line.  In  the  current  work  a  divergence  free  basis  is  constructed  for  a  standard  element  type 
and  curved  edges  are  allowed. 

Many  numerical  algorithms  have  been  developed  to  solve  this  problem.  Even  though  the  null  basis  is  not 
unique,  we  are  interested  in  those  which  are  sparse  and  banded.  This  means  small  support  basis  functions  for 
the  divergence  free  subspace  and  a  banded  sparse  transformation  matrix  for  the  dual  variable  method. 

The  turn  back  algorithm  is  a  method  for  computing  a  banded  and  sparse  null  basis.  It  was  proposed  by 
Topcu  [22]  to  compute  a  null  basis  with  a  profile  structure  for  equilibrium  matrices  in  structural  analysis. 
Kaneko,  Lawo  and  Thierauf  [23]  interpreted  this  algorithm  from  a  matrix  factorization  point  of  view.  Berry,  et 
al  [24]  refined  this  algorithm,  implemented  it  using  profile  data  structures  and  tested  it  on  several  structural 
problems.  Berry  and  Plemmons  [25]  have  implemented  this  algorithm  on  a  HEP  multiprocessor. 

Coleman  and  Pothen  [26]  have  designed  two  methods  to  obtain  a  sparse  null  basis:  the  first  one  computes 
a  fundamental  basis  (one  with  an  embedded  identity  matrix);  the  other  one  computes  a  triangular  basis  (one  with 
an  upper  triangular  matrix).  Both  algorithms  have  two  phases.  In  the  first  combinatorial  phase,  a  minimum 
dependent  set  of  columns  is  identified  by  finding  a  match  in  the  bipartite  graph  of  the  matrix.  In  the  second, 
numerical  phase,  nonzero  coefficients  in  the  null  vector  are  computed  from  this  dependent  set  Finding  the 
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sparsest  basis  for  the  null  space  of  an  underdetermined  matrix  was  shown  to  be  NP-hard  by  Coleman  and  Pothen 

[21]. 


We  have  developed  an  algorithm  which  yields  a  basis  for  the  divergence  free  subspace  associated  with  the 
standard  isoparametric  8-node  velocity  approximation  and  4-node  isoparametric  pressure  element.  This  basis  has 
minimal  support  (no  greater  than  9  elements),  can  be  efficiently  constructed,  and  curved  elements  are  allowed. 
The  relation  between  this  divergence  free  subspace  approach  and  the  dual  variable  method  has  been  established. 
The  turnback  algorithm  is  used  on  a  patch  of  elements.  This  work  is  contained  in  a  forthcoming  Ph.D.  thesis  of 
X.  Ye,  a  Ph.D.  student  of  Professor  Charles  Hall. 
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5.  Numerical  Solution  of  Transport,  Slightly  Diffusive  Transport  and  Fluid  Flow  Problems. 

A  central  theme  of  this  project  has  been  the  numerical  analysis  of  multidimensional  transport  and  slightly 
diffusive  transport  problems.  The  research  c.i  b^.  broadly  divided  into  the  following  (closely  interrelated)  areas: 

a.  Classification,  comparison  and  parameter  selection  for  difference 
schemes  for  2-D  convection-diffusion  problems. 

b.  Derivation  of  monotone  type  discretizations  for  boundary  value  problems. 

c.  Numerical  analysis  of  defect  or  deferred  correction  methods  for 
multi-dimensional  convection  diffusion  equations. 

d.  Numerical  and  analytical  studies  of  natural  convection  problems. 


a.  Classification  of  schemes  for  2-D  convection  diffusion  problems 
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In  Layton  [1],  structure  has  been  given  to  the  many  schemes  proposed  for  2-D  convection-diffusion  prob¬ 
lems  via  the  theory  of  modified  equations.  One  desirable  feature  in  a  numerical  methods  for  convection-diffusion 
problems  is  that  artificial  diffusion  operate  primarily  in  the  streamline  direction.  In  the  preliminary  study  [2]  it 
was  pointed  out  that  with  positive  type  stencils  for  general  (a,  b)  the  amount  of  "crosswind"  smearing  due  to 
artificial  viscosity  must  be  substantial  in  a  quantifiable  manner.  This  idea  was  more  fully  developed  (for  general 
stencils)  in  Layton  [1]. 

As  a  concrete  example  of  the  generally  applicable  theory  of  Layton  [1],  consider  the  celebrated  "skew- 
upwind"  scheme  (proposed  by  Raithby  in  1976)  which  was  a  precursor  of  the  currently  popular  streamline 
diffusion  methods.  In  [1]  it  was  shown  that  the  principal  axes  of  diffusion  of  this  scheme  are  not  aligned  with 
the  streamline  direction  unless  (a,  b)  is  aligned  with  the  mesh  [a  =  0,  b  =  0  or  a  =  b].  Futher,  in  [1,  example 
4]  it  is  shown  how  the  skew  upwind  scheme  should  be  "realigned".  The  resulting  scheme  is,  in  fact,  of  positive 
type  while  the  skew-upwind  scheme  is  not!  Numerical  experiments  in  [1,  section  4]  confirm  the  superiority  of 
this  new  scheme  over  both  the  original  skew-upwind  method  and  other  positive-type  upwinding  formulas. 

b.  Derivation  of  Monotone  Type  Discretizations  for  Boundary  Value  Problems 

The  well-known  stability  barrier  of  van  Veldhuizen  states  that  for  general  velocity  fields  (a,  b)  a  positive 
type  approximation  of  Lt  in  (1)  is  limited  to  first  order  accuracy.  In  Layton  [1]  it  was  shown  that,  even  among 
first  order  schemes,  for  general  (a ,  b )  positive  type  stencils  must  possess  substantial  cross-wind  diffusion  in  a 
quantifiable  sense.  Thus,  there  is  also  a  barrier  in  the  spurious  cross-wind  diffusion  in  these  schemes.  Further, 
for  second  order  elliptic  equations  with  cross-derivative  (u^)  term  there  is  also  a  stability  barrier  (due  to  Green¬ 
span  and  Jain)  on  the  size  of  the  u ^  coefficient  for  the  existence  of  positive  type  discretizations,  (see  Layton 
[14]  for  a  more  complete  discussion). 

Thus,  the  sign  pattern  required  of  positive  type  stencils  is  clearly  loo  restrictive  for  accuracy.  On  the 
other  hand,  it  is  very  desirable  to  have  discretizations  which  satisfy  a  discrete  maximum  principle/preserve  posi¬ 
tivity.  Thus,  schemes  which  lead  to  monotone  type  (or  inverse  positive)  discretization  matrices  must  be  con¬ 
sidered. 

This  condition  on  the  stencil  is  global.  One  contribution  of  this  research  is  that  it  has  been  shown  how  to 
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correctly  "microlocalize"  this  condition. 

A  general  methodology  for  the  microlocal  construction,  which  is  based  upon  the  theory  of  regular  split¬ 
tings  was  derived  and  analyzed  in  Layton  [15].  To  describe  it,  we  shall  focus  on  the  convection-diffusion  prob¬ 
lem  (1).  Consider  first  the  case  e  =  0,  (we  will  return  to  the  0  <  e  "small"  case  soon).  One  basic  problem  is 
going  from  the  microlocal  problem  (pre-monotone  stencils)  to  the  global  problem  (variable  coefficient  prob¬ 
lems).  Loosely  speaking,  this  difficulty  arises  because  matrix  multiplication  is  row  by  column,  thus,  a  point-by¬ 
point  splitting  of  L*  into  "generalized  positive  type  decompositions”  will  not  induce  a  matrix  splitting  when 
boundaries  or  variable  coefficients  are  present  Two  principles  of  "microlocalization"  have  been  proposed  which 
resolve  this  difficulty.  The  first  method,  in  Layton  [15],  is  based  upon  a  novel  theory  of  difference  schemes, 
intermediate  between  positive  type  and  general  monotone  type  schemes.  Positive  Averaged  Stencils  (PAS). 
These  were  introduced  in  Layton  [14]  where  the  relevant  theory  is  developed. 

Thus,  positive  averaged  stencils  are  one  way  to  pass  from  local  constructions  of  Lh  to  global  conditions 
on  (L*)-1!  Interestingly,  all  of  the  classical  monotone  type  discretizations  of  second  order  problems  are  positive 
averaged  stencils,  with  one  exception  in  Layton  and  Morley  [17].  The  extra  degrees  of  freedom  in  PAS  and 
relaxation  of  the  sign  pattern  can  be  used,  for  example,  to  increase  the  accuracy  in  methods. 

The  treatment  of  boundary  conditions  inside  the  P.A.S.  framework  is  not  difficult  To  preserve  monotoni¬ 
city  in  the  non-periodic  case  we  use,  at  the  points  adjacent  to  the  boundary,  a  lower  order  discretization.  This 
discretization  is  determined  by  computing  the  operator  splitting  Pk  Ah  restricted  to  those  points.  In  this  fashion 
we  generate  one  "compatible"  with  the  interior  equations. 

To  analyze  the  singularly  perturbed  case  a  "collage"  type  result  was  proven  in  Layton  [15]  monotone  (or 
inverse  positive)  matrices.  The  "collage  lemma"  of  Layton  [15]  gives  a  result  which  has  a  degree  of  uniformity 
in  the  dimension  of  the  matrices  involved. 

The  second  principle  of  microlocalization  was  presented  in  Layton  [14].  It  involves  differencing  variable 
coefficient  problems,  not  via  the  frozen  coefficient  problems  at  each  point,  but  rather  selectively  values  of  the 
coefficients  at  the  meshpoint  as  well  as  its  nearest  neighbors.  This  procedure  is  quite  general  and  systematic. 
The  derivation  of  the  variable  coefficient  stencil  can  be  summarized  as  follows: 
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Suppose,  in  the  constant  coefficient  case  a  stencil  Lh  possesses  a  "generalized  positive  type  decomposi¬ 
tion"  (in  the  sense  introduced  by  Brandi),  specifically 

[z/  =  A+A"  -  R  ,  R  >0  coef  ficientwise  , 

|A+  ,  A"  are  positive  type  difference  stencils  . 

Then,  the  correct  extension  of  Lh  to  variable  coefficient  problems  is  in  terms  of  the  frozen  coefficient  problems 
in  A+  and  A-  rather  than  in  Lk ! 

In  summary,  progress  on  monotone  type  discretizations  has  been  substantial.  The  question  of  derivation  of 
stencils  of  the  required  type  has  been  substantially  solved  and  promising  work  is  underway  on  the  remaining 
issue  (2-D  convection-diffusion).  The  'imuo-iocjli/ation”  question  is  completely  resolved  in  its  full  generality. 
Treatment  of  points  adjacent  to  computational  boundaries  has  been  successfully  addressed  in  several  examples 
by  ad-hoc  methods.  Additive  combination  of  monotone  type  stencils  with  positive  type  stencils  has  been 
resolved  in  the  case  of  a  "small"  perturbation  in  the  collage  lemma. 

c.  Numerical  Analysis  of  Deferred  Correction  Methods 

A  promising  algorithm  for  multidimensional  convection  problems,  originally  due  to  Hemker  [6],  based  on 
deferred  corrections  was  analyzed  in  Ervin  and  Layton  [3]  and  Axelsson  and  Layton  [4],  [5]  and  [18].  Addi¬ 
tional  work  included  extensive  numerical  experiments,  reported  in  Ervin  and  Layton  [13]  and  research  code, 
documented  in  Ervin  and  Layton  [16].  In  spite  of  intensive  computational  studies  of  (DCM),  rigorous  error 
analysis  was  lacking  for  problems  with  realistic  boundary  layers. 

In  Ervin  and  Layton  [3],  a  finite  difference  implementation  of  (DCM)  was  analyzed  for  the  two-point 
B.V.P.  arising  as  the  1-D  version  of  (1).  In  this  case,  uniform  in  e  (singular  perturbation  parameter)  conver¬ 
gence  of  optimal  orders  on  subdomains  was  proven.  Precise  estimates  near  layers  were  given  in  [3]  also. 

These  estimates  were  also  proven  to  be  sharp  in  numerical  experiments  reported  in  Ervin  and  Layton  [3]. 
The  local  error  analysis  via  cutoff  functions  was  carried  out  for  the  2-D  problem  in  Axelsson  and  Layton  [4]. 
This  analysis  focused  on  a  finite  element  interpretation  of  (DCM).  To  summarize  the  basic  result,  let  ft,  be  a 
nested  sequence  of  subdomains. 


fl  D  D  ^  •  =>  Ctj  , 


(3) 
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such  that  each  Q ,  does  not  admit  upstream  cutoff,  as  pictured  Figure  1 .  Further,  suppose 

(i)  the  characteristic  portion  of  DQ,  is  0(  .7T  In  (l/h )) 
distance  from 

(ii)  the  outflow  portion  of  dClj  is  0(/i  ln(l/h)) 
distance  from  j, 

(see  [4]  for  a  precise  formulation),  where  h  is  the  maximum  diameter  of  the  finite  element  triangulation  of  £}. 
The  basic  result  in  Axelsson  and  Layton  [4]  states  that  when  the  finite  element  space  approximates  smooth  func¬ 
tions  to  0 (hk)  in  //'(Q),  the  error  u  -  UJ  is  of  optimal  order,  uniformly  in  e  on  subdomains,  modulo  a  term 
which  is  of  infinite  order  accuracy,  nonuniformly  in  e. 

Theorem  Let  u  :=  true  solution  of  (1).  U1  :=  j‘h  finite  element,  defect  correction  approximation,  then  for 
any  s  >  0, 

i«  -  V(<v  -  -ty  +hk)  +  c2(£)e<5 . 

d.  Numerical  and  Analytical  Studies  of  Natural  Convection  Problems 

Convection-diffusion  equations  are  useful  models  of  more  complex  physical  processes,  such  as  coupled 
Navier-Stokes  system  in  free  or  natural  convection  problems,  studied  in  Boland  and  Layton  [7],  [8],  [12]  and 
Ermentrout,  Boland,  Hall,  Layton  and  Melhem  [9].  In  this  work,  a  finite  element  model  of  natural  convection  in 
an  enclosure  with  thermally  conducting  walls  is  studied,  based  upon  a  Boussinesq  model. 

In  a  preliminary  study,  the  stability  and  regularity  properties  of  the  Boussinesq  model  (4)  was  anal) zed. 
This  study  includes  the  discontinuous  coefficient  k  and  the  coupling  between  the  energy  and  the  Navier-Stokes 
equations.  This  preliminary  study  elucidates  expected  regularity  properties  of  the  solution  (u,  p,  T)  of  (4).  With 
this  knowledge  in  hand,  in  Boland  and  Layton  [7],  the  finite  element  method  for  the  time  dependent  coupled 
system  was  analyzed.  Error  estimates  were  proven  for  the  time  dependent  problem.  These  estimates  delineate 
both  expected  convergence  rates  under  minimal  data  restrictions  as  well  dependence  of  the  error  upon  the 
Prandtl  and  Rayleigh  (Pr ,  Ra )  numbers. 

In  a  companion  analysis  [12]  (incorporating  the  report  [8]),  we  consider  the  steady  state  problem.  Error 
estimates,  under  global  and  local  uniqueness  conditions,  are  proven  for  the  finite  element  model.  The  assump¬ 
tions  on  the  continuous  problem  are  also  proven  in  a  separate  section  of  [12].  One  interesting  auxiliary  result 
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proven  in  Boland  and  Layton  [19]  is  as  follows.  Consider  the  steady  state  problem  arising  before  subtracting  of 
the  boundary  conditions  specifically  y  =  0  with  inhomogeneous  Dirichlct  boundary  conditions  on  32Q  in  (4).  It 
was  shown  in  [12]  that  the  continuous  problem  always  has  at  least  one  steady  state  solution,  which  is  unique  for 
small  data  and  genetically,  locally  unique  for  large  data.  The  corresponding  question  can  be  asked  for  the  finite 
element  model:  does  the  finite  element  model  admit  steady  state  solution? 

Case  1:  dist(8Cle,  9Q)  >  0.  The  answer  here  is  yes. 


Case  2:  Linear  or  bilinear  elements  plus  an  angle  or  aspect  ratio  condition.  The  answer  here  is  yes,  via 
a  discrete  maximum  principle  argument  in  the  energy  equation. 

Case  3:  Quadratic  or  higher  elements  for  the  temperature.  In  this  case  the  following  condition  is  needed 
to  ensure  a  yes  answer:  cither  dist  (90,  90, )  >  0  or  O  =  O,  and:  the  local  mesh  density  near  the  vertical  sides 
of  90,  must  be  sufficiently  small,  in  a  quantifiable  sense,  depending  on  Pr  and  Ra  :  h]oc  =  C  Ra~]l4  near  90,. 
This  result  is  in  accordance  with  the  engineering  practice  of  grading  the  mesh  near  90,  as  in  (5).  In  fact,  a  for¬ 
mal  asymptotic  analysis  of  layers  in  (9)  suggests  that  the  boundary  layers  are  in  fact  within  0(/?a~1/4)  of  90, . 

This  idea  was  developed  in  a  different  direction  in  [9],  It  is  well  known  that  discretizations  of  nonlinear 
problems  frequently  introduce  spurious  solutions.  In,  e.g.,  Stephens  and  Shubin  [10],  it  was  shown  how  these 
spurious  steady  state  solutions  are  introduced  for  Burger’s  equation,  by  spatial  discretization.  Motivated  by  non¬ 
physical  dynamics  observed  in  the  finite  element  model  of  hollow  glass  blocks,  in  [9],  we  study  how  these  non¬ 
physical  steady  states  influence  the  dynamics  of  the  time  dependent  Burger’s  equation.  In  some  cases,  we  have 
rigorously  proven  that  a  Hopf  bifurcation  does  occur  in  the  discrete  system.  This  is  consistent  with  the  compu¬ 
tational  results  reported  in  Melhem  [11], 
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6.  Equilibrium  Manifolds  in  Computational  Mechanics 

Throughout  the  past  several  years,  one  aspect  of  our  research  program  concerned  the  computational  solu¬ 
tion  of  equilibrium  problems  as  they  arise  naturally  in  continuum  and  fluid  mechanics.  The  mathematical  models 
for  describing  these  problems  are  formed  by  nonlinear  equations,  including  algebraic,  differential,  or  integral 
equations,  all  of  which  involve  several  parameters,  and  hence  have  the  generic  form  F(z,X)  =  0  .  Here  z  varies 
in  some  space  Z  and  characterizes  the  state  of  the  physical  system  while  X  denotes  the  vector  of  parameters 
from  some  finite  dimensional  space  A. 

Our  research  in  this  area  has  been  based  on  a  consistent  use  of  the  fact  that  the  solution  set  of  such 
parametrized  equations  constitute  in  general,  a  differentiable  manifold  in  the  product  space  X  =  Z  xA.  Our 
results  so  far  certainly  show  the  value  and  power  of  this  geometric  approach. 

Our  work  concerns  essentially  the  following  topics: 

(1)  Development  of  methods  for  analyzing  various  features  of  equilibrium  manifolds  on  the  basis  of  a  com¬ 
puted  simplicial  approximation. 

(2)  Methods  for  the  computation  of  sub-manifolds  of  foldpoints  on  equilibrium  manifold. 

(3)  Study  of  estimation  techniques  for  the  errors  arising  in  the  various  computational  procedures  under  (1)  and 

(2). 

Some  details  of  the  current  work  are  given  in  the  following  sub-sections. 

i.  Methods  for  Analyzing  Features  of  Equilibrium  Manifolds. 

A  new  method  for  computing  simplicial  approximation  of  the  solution  manifolds  of  a  parametrized  equa¬ 
tion  has  been  developed.  An  essential  par!  of  the  method  is  a  constructive  algorithm  for  computing  moving 
frames;  that  is,  of  orthonormal  bases  of  the  tangent  spaces  that  vary  smoothly  with  their  points  of  contact  The 
approximation  process  uses  these  bases  together  with  a  chord  form  of  the  Gauss-Newton  method  as  a  corrector. 


to  compute  the  desired  nodes  of  the  simplicial  mesh  on  the  manifold. 
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The  process  has  been  implemented  for  two-dimensional  as  well  as  three-dimensional  manifolds  in  spaces 
of  arbitrary  dimension.  All  experience  with  these  codes  so  far  has  been  excellent.  In  particular,  it  turns  out  that 
the  method  is  extremely  efficient  and  in  fact  has  a  better  computational  complexity  than  most  continuation  pro¬ 
cedures.  For  example,  for  a  two-dimensional  manifold  a  typical  computation  produces  a  triangulation  involving 
1 14  triangles  with  19  evaluations  and  decompositions  of  the  Jacobian  of  the  mapping.  Work  is  continuing  on  the 
process.  In  particular,  we  are  currently  experimenting  with  adaptive  approaches  for  handling  situations  involving 
large  changes  of  the  curvature  of  the  manifold. 

Any  simplicial  approximation  computed  with  this  method  provides  a  wealth  of  information  about  the 
equilibrium  manifold.  This  opens  up  numerous  possibilities  for  determining  various  features  of  such  manifolds. 

A  first  possibility  is  of  course  to  provide  a  graphical  representation  of  parts  of  the  manifold.  Such  an 
opportunity  has  not  existed  before  and  already  simple  examples  indicate  what  new  insight  can  be  obtained  in 
this  way.  At  the  same  time,  effective  graphics  techniques  for  visualizing  a  p-dimensional  manifold  defined  in  a 
high  dimensional  space  are  by  no  means  readily  available.  During  the  past  year  a  set  of  post-processing  routines 
have  been  developed  which  generate  a  sequence  of  slices  of  a  simplicial  approximation  of  a  three-dimensional 
manifold.  Suitable  projections  of  these  two-dimensional  slices  can  then  be  represented  graphically  by  means  of 
standard  programs,  such  as  MOVIE.BYU.  By  viewing  these  slices  in  sequence  one  has  the  impression  of 
"flying"  through  the  3-dimensional  muiuloiu  anu  tins  provides  often  some  unique  insight  into  its  features.  Work 
is  now  continuing  to  incorporate  plots  of  prescribed  contour  surfaces  into  these  graphical  representations. 

Graphical  representations  are,  of  course,  only  one  way  to  analyze  the  various  features  of  our  manifolds. 
The  wealth  of  computed  data  offer  numerous  possibilities  for  the  application  of  post-processing  routines  to  cal- 
ulate  approximations  of  various  specific  features. 

A  particularly  important  feature  is  the  curvature  tensor  of  the  manifold  M  which  has  many  important 
applications.  The  numerical  calculation  of  the  curvature  tensor  from  the  data  made  available  by  the  moving 
frame  algorithm  was  set  as  one  of  the  goals  of  the  previous  proposal.  The  differential  geometric  and  numerical 
difficulties,  the  latter  being  essentially  due  to  the  fact  that  curvature  depends  upon  second-order  quantities,  have 
been  successfully  overcome  in  [1].  The  corresponding  algorithm  has  provided  surprisingly  good  results.  More¬ 
over,  an  unexpected  by-product  of  this  method  was  found  in  the  form  of  the  apparently  first  efficient  and 
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rigorous  technique  for  the  calculation  of  all  the  bifurcated  branches  of  solutions  in  bifurcation  problems  involv¬ 
ing  a  higher  dimensional  null-space.  Up  to  now,  only  the  one-dimensionai  null-space  case  had  been  considered 
(and  so,  extensively)  in  the  numerical  literature. 

As  noted  earlier  the  availability  of  some  approximation  of  the  curvature  tensor  has  many  important  appli¬ 
cations.  One  of  them  concerns  the  question  of  estimating  the  distance  from  a  given  point  x  on  M  to  some  fold- 
point  x*  of  that  manifold.  As  part  of  the  methods  for  computing  foldpoints  discussed  in  the  next  subsection  we 
have  developed  a  method  for  obtaining  such  estimates.  In  essence,  a  foldpoint  on  M  is  defined  by  the  property 
that  one  or  several  normal  vectors  of  M  at  that  point  are  parallel  to  the  given  natural  parameter  space  A.  This 
suggests  the  computation  of  the  smallest  principal  angle  between  the  (n-p)-dimensional  normal  space  rgeDF(x)* 
and  the  p-dimensional  natural  parameter  space  A.  These  principal  angles  arc  readily  computed,  but,  of  course, 
without  any  further  information  they  do  not  provide  the  desired  estimate  of  llx-x*ll.  Here  the  approximate  curva¬ 
ture  tensor  can  be  used  very  effectively.  This  provides,  for  the  first  time,  a  computable  estimate  for  this  impor¬ 
tant  distance.  Results  so  far  with  the  technique  arc  very  encouraging  but  further  work  is  needed  to  make  the 
method  generally  applicable. 

The  computed  simplicial  approximations  also  provide  an  effective  tool  for  deterr.  ning  the  approximate 
location  of  the  foldset  of  the  manifold  M.  For  this  one  has  to  monitor  the  orientation  of  the  projection  of  the 
tangent  basis  at  the  nodal  points  onto  the  parameter  space  A.  If  there  is  a  change  in  this  orientation  then  we  are 
in  the  neighborhood  of  a  foldpoint;  but  the  converse  is  not  necessarily  true,  since  not  every  foldpoint  can  be 
detected  this  way.  The  orientation  can  be  characterized,  for  example,  by  the  determinant  of  the  projected  basis  in 
the  parameter  space.  Thus,  in  essence,  if  we  plot  lines  of  constant  determinant  values,  then  lines  of  zero  deter¬ 
minant  values  are  approximations  of  the  desired  foldlines.  Of  course,  determinants  are  not  always  the  best  way 
to  proceed  here,  and  we  are  currently  studying  alternate  techniques  for  the  effective  determination  of  the  desired 
orientation  changes  and  for  combining  the  technique  with  the  above  mentioned  method  for  estimating  the  dis¬ 
tance  to  the  desired  foldpoint(s). 

ii.  Methods  for  Computing  Foldpoints. 

One  of  the  central  computational  problems  in  connection  with  the  equilibrium  manifold  M  of  a  given 


parametrized  equation  is  the  determination  of  the  foldsets  of  M.  These  sets  are  of  special  importance  in  stability 
considerations  for  the  equilibrium  problems  under  study. 

As  mentioned  in  the  previous  sub-section,  the  availability  of  the  simplicial  meshes  opens  up  a  unique  way 
of  detecting  the  presence  of  such  foldpoints.  The  points  on  the  simplicial  approximation  also  provide  approxima¬ 
tions  for  these  foldpoints  from  where  local  iterative  processes  can  be  started  for  their  accurate  calculation.  All 
such  iterative  methods  require  a  suitably  augmented  system  of  equations  for  which  the  Jacobian  is  non-singular 
at  the  desired  foldpomt.  For  some  time  now,  the  study  of  various  classes  of  such  augmenting  schemes  has  been 
an  active  topic  of  research  under  this  project  In  particular,  in  [2]  we  consider  augmentations  of  the  form 

F(x)  =  Q 

cTUj  -  0  j= 1,  •  •  ,  p 

where  the  u,  form  an  orthonormal  basis  of  the  tangent  space  at  the  point  x  and  c  is  a  vector  in  the  intersection 
of  the  normal  space  at  that  point  with  the  parameter  space  A.  This  augmentation  was  shown  to  be  feasible  at 
singular  points  of  socallcd  type  <p,l,l)  which  includes  the  standard  turning  points,  simple  bifurcation  points,  and 
also  various  other  types  of  singularities. 

The  crucial  point  in  the  utilization  of  this  augmentation  is  the  choice  of  the  ’’unfolding’’  vector  c.  During 
the  past  year  a  new  method  has  been  developed  to  obtain  such  a  vector  c  when  x  is  an  approximation  of  some 
foldpoint  x*  of  the  manifold  [3],  The  technique  is  closely  related  to  the  approach  mentioned  in  the  previous 
sub-section  for  estimating  the  distance  I  Ix-x*  I  I.  A  principal  idea  is  the  use  of  the  so-called  gap  between  the 
tangent  space  at  x  and  the  parameter  space  A.  As  part  of  this  approach  it  has  become  possible  to  generate  a  new 
augmentation 

F(x)  =  0 

g,(x)  =  0  ,  j= I,  ■  •  ,  p 

with  the  desirable  property  that  the  finite  difference  approximations  of  the  derivatives  of  the  gj  can  be  obtain 
directly  from  the  Jacobian  of  F.  Moreover,  for  paths  of  foldpoints  this  particular  choice  of  the  augmentation 
allows  us  to  determine  the  approximate  direction  of  a  path  of  foldpoints  on  the  manifold.  With  this  we  have 
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obtained  a  new  method  for  continuing  along  such  paths  of  foldpoints  which  turns  out  to  be  considerably  more 
efficient  than  other  methods  known  so  far.  More  generally,  the  new  augmentation  has  also  opened  up  the  possi¬ 
bility  of  determining  the  tangent  directions  of  a  higher  dimensional  sub-manifold  of  foldpoints.  This,  in  turn, 
provides  for  the  application  of  the  simplicial  approximation  algorithm  to  such  sub-manifolds.  Such  "triangula- 
tions"  of  foldsets  have  never  been  obtained  before  and  they  offer  promising  new  approaches  for  the  analysis  of 
stability  problems. 

iii.  Error  Estimations 

There  are  numerous  error  questions  arising  in  the  different  processes  for  computing  foldpoints  on  a  solu¬ 
tion  manifold.  During  the  past  year  we  have  begun  a  systematic  study  of  the  sources  of  these  errors  and  of 
methods  for  estimating  them.  This  work  has  addressed  the  following  computational  tasks  related  to  the  determi¬ 
nation  of  foldpoints. 

(a)  Detection  problem: 

This  important  problem  may  be  phrased  a  ;  follows:  Given  a  set  (x1 . xk )  of  points  and  appropriate  addi¬ 

tional  information  about  the  mapping  at  these  points,  determine  whether  some  foldpoint  of  the  manifold  M  is  in 
a  neighborhood  of  this  set. 

This  is  in  essence  a  question  of  determining  how  far  a  given  point  is  from  a  solution  of  a  suitable  aug¬ 
mented  system  of  equations.  Evidently,  on  the  basis  of  information  at  finitely  many  points  alone  we  cannot 
expect  to  obtain  a  guaranteed  inclusion  result;  at  least  some  information  about  the  mapping  in  a  suitable  neigh¬ 
borhood  of  the  given  set  of  points  is  needed.  This  is  already  seen  lor  very  simple  cases.  Suppose  that 
\=span (v )  and  that  x]  and  x2  arc  two  successive  points  produced  by  a  continuation  procedure.  If  both  points 
belong  to  the  same  connected  component  of  the  continuation  path  and  u‘  and  u2  arc  the  corresponding 
(oriented)  tangent  vectors  then  sgn(vT ul)*sgn(vT u2)  implies  that  there  is  a  limit  point  of  odd  order  between 
these  poirts.  Evidently  without  further  information  about  the  mapping  we  cannot  verify  that  both  points  belong 
to  the  same  connected  component  of  the  path  and  without  that  assumption  the  statement  need  not  hold.  On  the 
other  hand,  suppose  that  we  only  know  that  the  two  points  are  on  the  path  and  that  we  have  oriented  the  tangent 
vectors  by  the  requirement  (u  l)T(u2)  >  0.  Now,  it  turns  out  that  when  a  certain  determinant  changes  sign 
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between  the  points  and  the  points  are  not  too  far  apart,  then  they  cannot  belong  to  the  same  connected  com¬ 
ponent  and  hence  there  must  be  a  bifurcation  point  between  them. 

These  observations  extend  to  the  general  case  of  a  p-dimensional  manifold.  Once  again  the  availability  of 
a  simplicial  approximation  of  the  manifold  and  of  some  derived  quantities,  as,  for  instance,  an  approximate  cur¬ 
vature  tensor,  provides  here  some  of  the  needed  information  for  the  development  of  methods  for  handling  the 
detection  problem. 

(b)  Approximation  error 

For  a  given  approximation  x  of  a  foldpoint  x*  of  M  estimate  the  distance  I  Ix-x*  I  1  between  these  points. 
This  depends  on  the  measure  of  the  distance  between  the  points;  at  the  same  time  some  information  about  the 
mapping  is  needed  to  derive  estimates  of  the  distance.  Here,  as  noted  in  subsection  1  above  the  size  of  the  smal¬ 
lest  principal  angle  between  the  m-dimensional  normal  space  of  the  manifold  at  the  point  x  and  the  p- 
dimensional  natural  parameter  space  A  gives  some  information  about  this  error.  These  principal  angles  are 
readily  computed,  but  again  information  about  the  mapping  is  needed  to  derive  from  this  an  actual  estimate  of 
I  lx-x*  I  I.  As  mentioned  earlier,  an  approximate  value  of  the  local  curvature  tensor  proves  very  useful  in  this 
connection. 

(c)  Foldpoint  computation 

We  have  discussed  already  methods  for  the  computation  of  specific  foldpoints  on  our  manifold.  In  all 
cases  a  suitably  augmented  system  of  equations  is  formed  which  then  allows  the  application  of  a  standard, 
locally  convergent  iterative  process.  The  behavior  of  this  process  is  controlled  by  the  error  of  the  initial  approxi¬ 
mation,  the  termination  error  of  the  method,  and  the  influence  of  round-off  on  the  reliability  of  the  result.  The 
initial  approximation  error  was  already  discussed  under  (a)  above.  But  the  application  of  the  iterative  process 
provides  some  additional  data  which  permit  a  Juscr  estimation  of  that  error.  In  connection  with  the  termination 
error  we  arc  not  only  interested  in  the  question  how  close  the  computed  point  is  to  the  manifold,  but  how  well  it 
actually  approximates  the  desired  foldpoint  Thus  we  are  again  back  to  the  previous  problem.  But  once  again, 
the  behavior  of  the  iterative  process  provides  some  information  that  was  not  available  before  and  that  appears  to 
help  in  gaining  a  closer  estimate. 
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1.  Differential  Algebraic  Equations  (DAE’S)  in  Mechanics 
Many  problems  in  mechanics  are  of  the  generic  form: 

F(x)  =  0  (1) 

A(x)x=G(x)  (2) 

where  x  =  (z , X),  z  is  from  a  state  space  Z  and  X  is  from  a  p  dimensional  parameter  space  A.  Equation  (1)  is 
a  set  of  nonlinear  algebraic  equations  and  (2)  is  a  set  of  differential  equations,  hence  the  system  (l)-(2)  is 
termed  a  DAE. 

An  alternate  interpretation  comes  from  the  observation  that  (1)  defines  a  manifold  M  of  dimensions  p 
in  the  combined  space  X  =  ZxA  and  that  (2)  is  a  differential  equation  on  M;  hence  the  system  (l)-(2)  is  also 
referred  to  as  a  differential  equation  on  a  manifold  (DEM). 

Two  applications  of  interest  to  the  investigators  in  which  DEM’s  or  DAE’s  occur  are: 

(i)  Incompressible  Fluid  Flow.  The  continuity  equation  is  an  algebraic  constraint  of  the  form  (1)  and  the 
time-dependent  Navier-Stokes  equation  is  a  differential  equation  of  the  form  (2). 

(ii)  Punch  Stetching  of  Sheet  Metal.  The  principle  of  virtual  work  provides  a  force  equilibrium  equation 
which  defines  an  equilibrium  manifold  upon  which  one  seeks  solutions  to  differential  constitutive  laws  of 
the  form  (2). 

i.  Elastic-Plastic  Deformation 
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In  [12],  a  new  numerical  method  was  presented  for  computer  simulations  of  punch  stretching  of  sheet 
metal.  Most  current  approaches  to  finite  element  modeling  of  large  deformation,  elastic-plastic  sheet  metal 
forming  use  a  rate  form  of  the  equilibrium  equations  and  then  must  correct  at  each  time  step  to  insure  that 
equilibrium  is  satisfied.  Such  methods  are  referred  to  as  incremental  methods  [3].  The  new  method,  a  DEM 
approach  discretizes  the  more  fundamental  equilibrium  equations  in  non-rate  form  and  insures  equilibrium  of 
forces  at  each  time  step.  Formulating  the  problem  as  a  DEM  or  DAE  also  allowed  for  solution  of  the  discre¬ 
tized  system  using  off-the-shelf  software  such  as  LSODI  [14].  Numerical  experimentation  indicated  that  the 
DEM  approach  was,  in  many  cases,  computationally  more  efficient  than  the  incremental  approach. 

Differential-algebraic  systems  (DAE’s)  are  usually  solved  by  means  of  software  developed  originally  for 
the  solution  of  initial  value  problems  of  ordinary  differential  equations  (ODE’s).  This  approach  has  become 
widely  accepted  ever  since  it  was  proposed  by  Gear  [4],  Such  codes  as  LSODI  [14]  and  DASSL  [15]  are  based 
on  this  approach  in  which  the  algebraic  equations  are  viewed  as  degenerate  ODE’s. 

An  alternative  approach  is  the  DEM  approach  mentioned  earlier,  (see  also  [5]).  The  DAE  is  viewed  as  a 
differential  equation  on  a  manifold;  the  latter  determined  by  the  algebraic  equations.  Hence  we  seek  a  curve 
passing  through  the  initial  point  whose  tenant  conforms  to  the  direction  field  defined  by  the  differential  equa¬ 
tions.  This  interpretation  allows  for  the  application  of  the  differential-  geometric  techniques  discussed  in  earlier 
sections.  In  particular,  it  suggests  to  new  methods  for  the  solution  of  DAE’s  especially  in  cases  where  the  ODE 
methods  no  longer  work  very  well. 

Our  experience  in  [12,13]  gave  strong  indications  that  the  DAE’s  arising  in  finite  element  discretizations 
of  clasto-viscoplastic  problems  can  not  always  be  solved  efficiently  by  ODE-soflware,  such  as  LSODI.  Problems 
appear  to  be  related  to  the  inability  of  the  predictor  to  handle  sudden  changes  in  external  forces  due  to  punch  or 
die  contact,  and  to  the  onset  of  plasticity. 

Computerized  mathematical  models  which  produce  accurate  and  efficient  simulations  of  sheet  metal  form¬ 
ing  processes  arc  highly  useful  for  example  in  the  automobile  industry  as  a  means  of  predicting  the  success  or 
failure  of  a  punch-die  design,  or  evaluating  formability  properties  of  new  materials.  The  differential  equations 
on  a  manifold  (DEM)  method,  coupled  with  state-of-the-art  software  packages  designed  to  integrate  differential 
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algebraic  equations  (DAE’s),  has  proven  to  be  such  a  model. 

The  basic  equations  governing  the  solution  of  large  deformation  plasticity  problems  in  continuum  mechan¬ 
ics  are  the  equilibrium  equations  and  the  constitutive  equations.  In  discretized  form,  the  equilibrium  conditions 
become  a  system  of  nonlinear  algebraic  cquauons,  while  the  elastic-plastic  constitutive  equations  are  a  system  of 
ordinary  differential  equations  (ODE’s).  Two  general  approaches  to  the  solution  of  such  differential-algebraic 
equations  (DAE’s)  are  typically  used  in  displacement-based  finite  element  methods.  In  one  approach,  the  time 
derivative  of  the  equilibrium  equations  is  calculated  and  the  DAE  becomes  an  ODE  in  nodal  displacements.  In 
another,  the  equilibrium  equation  is  left  unchanged  and  during  each  time  iteration  of  the  nodal  displacements, 
the  constitutive  equations  are  integrated  to  yield  the  stress,  which  is  then  used  to  check  for  equilibrium. 

The  Differential  Equation  on  a  Manifold  (DEM)  method  as  presented  in  [12]  is  similar  to  the  latter 
approach  in  that  the  original  equilibrium  equations  are  used,  however,  as  in  the  mixed  methods  the  stresses  as 
well  as  the  displacements  have  finite  element  approximations.  The  equilibrium  equations  are  viewed  as  a  mani¬ 
fold  in  the  solution  space  on  which  the  solutions  to  the  ODE  (constitutive  equations)  must  be  traced.  State-of- 
the-art  numerical  software  is  then  used  to  solve  the  resulting  DAE.  Equilibrium  is  satisfied  at  each  time  step. 

The  DEM  method  was  applied  successfully  to  axisymmetric  hydrostatic  bulging  [12]  and  to  axisymmetric 
and  plane  strain  sheet  metal  forming  problems  [13].  As  compared  to  previous  results  [3],  the  DEM  method  pro¬ 
duced  accurate  solutions  for  rate-sensitive  materials  for  computational  running  times  which  range  from  6  to  26 
times  less  than  for  previous  methods.  The  software  package  LSODI  [14]  was  used  for  the  numerical  integration 
of  the  DAE’s.  These  implementations  of  'h:  DEM  method  for  punch  stretching  resulted  in  an  efficient  and 
accurate,  but  somewhat  non-robust  algorithm. 

In  order  to  assess  better  the  expected  behavior  of  the  solutions,  a  study  of  the  existence  and  uniqueness  of 
nonlinear  DAEs  was  begun.  In  [21]  such  a  theory  was  presented  for  first  and  second  order  systems  of  the  form 

F  ,00  =  0 

F2(x,  x',z)  =  0 

and 

F  ,00  =  0 

F2(x,  x' ,  x",  z)  =  0  , 
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respectively.  For  these  systems  it  was  shown  that  they  induce  smooth  vector  fields  on  certain  open  subsets  of 
the  space  of  their  differentiable  variables  r  Moreover,  in  the  first  order  case  these  vector  fields  are  tangential  to 
the  constraint  manifolds  (x;  F  fx)  =  constant}  while  for  the  second  order  equations  they  are  tangential  to  the 
tangent  bundles  of  these  manifolds.  This  in  turn  provides  directly  for  the  desired  existence  and  uniqueness  of 
these  DAEs.  The  theory  has  also  opened  up  the  development  of  a  new  local  parametrization  approach  for  the 
computational  solution  of  these  systems.  This  generalizes  earlier  developed  techniques  for  the  Euler-Lagrange 
equations  arising  in  constrained  multibody  dynamics. 

ii.  Bifurcation  Phenomena 

Our  computational  experience  with  sheet-metal  forming  has  also  given  strong  indication  that  in  the  solu¬ 
tion  of  DAE’s  there  are  situations  in  which  singularities,  such  as  bifurcations,  occur.  In  particular,  certain  non¬ 
linear  DAE’s  may  exhibit  multiple  solutions  in  some  domains.  Two  examples  of  such  a  situation  are  the  follow¬ 
ing: 

(a)  A  geometrically  symmetric  plane  strain  punch  problem,  such  as  we  considered  in  [12,  13],  was  modeled 
successfully  using  a  natural  symmetry.  However,  when  this  symmetry  was  ignored  the  solver  (in  this  case 
LSODI)  produced,  what  appeared  to  be,  two  distinct  solution  curves  after  the  onset  of  plasticity  depending 
on  the  particular  problem  parameters.  One  of  these  curves  corresponded  to  the  expected  symmetric  solu¬ 
tion  while  the  other  one  was  clearly  un-symmctric,  which  suggests  the  occurrence  of  a  symmetry  breaking 
bifurcation  point. 

(b)  An  extension  of  the  DEM  approach  to  in-plane  stretching  by  Punch  [6]  led  to  a  similar  situation  in  which 
LSODI  failed  to  follow  a  geometrically  symmetric  solution  curve  in  favor  of  one  which  was  non- 
symmctric. 

These  examples  point  to  the  necessity  of  a  closer  study  of  the  singularities  of  DAE’s.  It  appears  that  up  to 
now  very  few  results  are  available  in  this  area. 

The  recent  theoretical  work  [10],  gives  n  complete  and  rigorous  analysis  of  the  dynamics  in  the  vicinity  of 
a  generic  singular  point.  The  results  in  [10]  explain  why  such  deviations  from  the  expected  symmetric  solution 
occur,  which  are  due  to  the  slightest  breaking  of  a  perfectly  symmetric  problem.  Numerically,  such  a  breaking 
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may  be  caused  e.g.  by  roundoff.  More  interestingly,  since  actual  problems  are  never  perfectly  symmetric,  these 
results  show  that  the  true  solution  may  well  eventually  lose  symmetry  completely,  because  the  symmetric 
configuration  becomes,  loosely  speaking,  "unstable"  beyond  the  singularity.  This  demonstrates  that  enforcing 
symmetry  of  the  solution  by  reducing  the  size  of  the  system,  although  tempting  numerically,  is  not  a  safe  pro¬ 
cedure  if  a  singularity  is  encountered,  lor  uien  uic  calculated  solution  is  no  longer  significant  beyond  the  singu¬ 
larity.  This  phenomenon  is  largely  (and  probably  totally)  unknown  and  its  discovery  could  have  considerable 
importance  in  problems  in  which  symmetry  of  the  solution  is  routinely  taken  for  granted  because  it  is  "obvi¬ 
ous.".  From  now  on,  consideration  should  be  given  to  the  nonobvious  and  unsuspected  fact  that  the  slightest 
breaking  of  symmetry  has  drastic  effects  on  the  solution  in  fixed  finite  time. 

In  the  case  of  the  punch  problem  described  in  [13],  a  careful  monitoring  of  the  eigenvalues  by  Y.  Huang, 
a  Ph.D.  student  of  Professors  Hall  and  Rabier,  has  shown  that,  indeed,  a  singularity  of  the  expected  type  is 
involved.  This  has  been  observed  using  various  mesh  sizes  for  the  discretization  and  only  little  dependence  on 
the  critical  time  has  been  observed.  A  summary  of  these  calculations  with  expanded  comments  is  currently  in 
preparation  ([11]). 

Of  course,  the  understanding  of  the  dynamics  in  the  vicinity  of  generic  ("standard"  in  the  terminology  of 
[10])  singular  points  of  differential  equations  allows  one  to  ask  precise  corresponding  numerical  questions,  with 
a  natural  priority  given  to  the  development  of  reliable  techniques  for  the  characterization  and  computation  of 
such  points.  One  difficulty  in  this  respect  is  that  the  approach  taken  in  [10]  is  essentially  theoretical  and  has  no 
numerical  counterpart.  However,  striking  analogies  between  these  standard  singular  points  and  the  much  more 
familiar  turning  points  in  transcendental  equations  make  it  plausible  that  appropriate  variants  of  methods  for  the 
calculation  of  turning  points  exist,  that  can  successfully  be  used  for  the  determination  of  singular  points. 

The  investigation  of  problems  exhibiting  symmetries  has  led,  indirectly  and  unexpectedly,  to  new  results 
in  problems  of  bifurcation  involving  symmetry  (see  [20])  and  having  consequences  in  abstract  group  theory  itself 
([21]).  These  results  also  have  important  potential  applications  to  the  numerical  treatment  of  such  problems 


which  are  currently  under  study. 
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8.  Bifurcation  Theory  and  Singular  ODE’s. 

In  connection  with  differential  algebraic  equations  (DAE’s),  a  theoretical  analysis  of  singular  differential 
equations  has  been  started  in  [5].  The  problem  is  to  determine  the  structure  of  the  set  of  solutions  to  an  initial 
value  problem 

|A(x)x  =  G(x)  , 

(*Oo)  =  x0  ,  ^ 

when  the  n  x  n  matrix  A(x)  is  singular  at  x  =  x0.  The  above  formulation  is  suitable  for  the  analysis  of  singu¬ 
larities  occurring  in  DAE’s  as  well,  a  framework  in  which  the  problem  has  received  a  great  deal  of  attention 
from  singular  perturbation  theory  specialists.  The  approach  taken  in  [5]  follows  ideas  of  bifurcation  theory  and 
has  permitted,  apparently  for  the  first  time,  to  eliminate  the  usual  but  restrictive  assumption  that  A  (x)  is  singular 
on  an  entire  neighborhood  of  x0. 

In  actual  applications,  equation  (1)  generally  does  not  exhibit  a  singularity  at  the  beginning  of  evolution. 
However,  it  is  quite  possible  that  singularities  appear  at  some  later  time  along  a  trajectory.  To  describe  what 
happens  next,  one  must  analyze  problems  of  the  form  (1)  with  singular  A(x0).  As  is  shown  in  [5],  existence  of 
a  singularity  along  a  trajectory  is  not  affected  by  small  perturbations  of  the  (nonsingular)  initial  condition:  it  is 
a  stable  and  hence  important  phenomenon. 

Our  recent  work  on  singular  differential  equations  has  followed  two  directions.  First,  we  have  derivel 
reliable  algorithms  for  the  calculation  of  the  singular  points  called  "standard"  in  [4].  These  points  are  of  pri¬ 
mary  interest  because,  according  to  [5],  most  singular  points  are  ot  this  type  in  the  generic  case. 

Starting  with  a  problem  of  the  form  (I)  with  nonsingular  A  (xa ),  we  assume  existence  of  an  unknown  time 
f.  >  0  such  that  A(x(f. ))  is  singular  and  propose  methods  for  calculating  t,  and  x.  =  x((>)  in  the  assumption 
that  x.  is  a  standard  singular  point 

The  need  for  such  methods  can  be  hinted  from  the  fact  that  the  simple  approach  consisting  in  writing  (1) 
in  the  explicit  form 

fx  =  A~l(x)G(x)  , 

[x(0)  =  xo 
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although  legitimate  for  t  €  [0,  /. ),  leads  to  a  problem  in  which  the  matrix  A(x(t ))  becomes  singular  as  t 
approaches  t. .  This,  of  course,  may  result  in  inaccuracies  in  the  numerical  calculation  of  both  r.  and 
x.  =  x  {t. ).  It  is  important  to  point  out  that  the  singularity  cannot  be  jumped  over.  Indeed,  the  analysis  in  [5] 
reveals  that  a  trajectory  encountering  a  standard  singular  point  must  necessarily  terminate  at  the  point:  the  solu¬ 
tion  x(t)  cannot  be  extended  continuously  for  t  >  t,.  As  far  as  numerical  aspects  are  concerned,  this  forces  the 
matrix  A(x(f))  to  remain  in  a  state  of  quasi -singularity  without  possible  escape  as  t  approaches  or  exceeds  /. 
(numerically,  t  may  exceed  /.  despite  the  fact  that  this  is  theoretically  impossible).  This  emphasizes  the  impor¬ 
tance  of  having  ad-hoc  algorithms  for  the  accurate  calculation  of  standard  singular  points. 

Such  algorithms  have  been  developed,  which  reduce  the  problem  to  solving  a  nonsingular  ODE  and  scalar 
algebraic  equations.  Part  of  the  procedure  uses  a  variant  of  a  minimal  augmentation  technique  introduced  by 
Griewank  and  Rcddien  [2]  and  extended  by  Rabier  and  Rcddien  [4],  in  connection  with  calculation  of  singulari¬ 
ties  in  bifurcation  problems.  The  possibilitv  of  modifying  methods  of  numerical  bifurcation  theory  to  handle 
singular  differential  equations  was  anticipated  in  the  previous  research  proposal.  A  related  report  [6]  is  currently 
being  prepared. 

The  second  component  of  the  work  recently  accomplished  in  the  domain  of  singular  differential  equations 
deals  with  their  relevance  in  real  life  problems.  Some  of  the  authors  of  this  proposal  had  been  suspecting  for  a 
while  that  numerical  inconsistencies  observed  in  a  problem  of  metal  punching  (see  [1])  were  due  to  the  presence 
of  a  singularity.  This  has  now  been  fully  corroborated  by  numerical  calculations  based  on  the  previously  men¬ 
tioned  algorithms  for  the  calculation  of  singular  points.  A  report  by  Hall,  Huang  and  Rabier  ([3])  on  this  prob¬ 
lem,  and  its  importance  and  implicauuiis,  in  picparation.  Figure  1  illustrates  the  behavior  of  X,  the  smallest 
eigenvalue  (in  modulus)  of  the  Jacobian  of  the  DAE  (differential  algebraic  equation)  resulting  from  the  DEM 
modelling  of  a  symmetric  punch  stretching  of  sheet  metal  problem  described  in  (1],  As  \  tends  to  zero  in  the 
vicinity  of  28  secs,  the  numerical  solution  loses  symmetry.  It  has  been  verified  numerically,  that  the  purely 
symmetric  problem  in  fact  has  a  Jacobian  which  is  singular  near  this  value  of  time.  (See  [3]  for  more  details.). 
For  various  finite  element  meshes  behavior  is  qualitatively  similar. 


Figure  1. 
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